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ABSTRACT 



Composite reliability is strongly dependent on the fiber strength 
distribution. Current methods of gathering statistics through single fiber failure 
methods are inefficient and costly. This thesis develops a testing method from 
which the fiber statistics in the form of Weibull parameters can be accurately 
extracted from bundle failure tests. 

The values obtained from the bundle experiment as compared to known 
single fiber test bench mark parameters were practically indistinguishable. The 
confidence of the results stems from a thorough analysis of the bundle 
mechanisms and minimized contaminates which can disturb the strength 
distribution. 
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I. INTRODUCTION 



A composite, by definition, is any two or more materials combined on a 
macroscopic scale to form a useful material [Ref 1]. Today, in one form or 
another, composites are being used on every level of our society. This thesis 
will focus on a class of composites known as "fibrous composites," today these 
are typically used in high tech applications, often high performance military 
aircraft. Fiber composites consist of very thin fibers with high aspect ratios and 
high strengths imbedded in a matrix, commonly a polymeric matrix such as 
epoxy. Glass/epoxy is an example of a frequently utilized composite mixture 
that has been in use for many years. 

One of the many advantages this class of material provides is high 
strength-to-weight, though due to a lack of thorough understanding in failure 
mechanics, composites have been limited to secondary load carrying structural 
designs such as skins and control surfaces. In recent years (last 20), composites 
are expanding into a much wider role. 

A. REQUIREMENTS IN ADVANCED STRUCTURES 

Since the introduction of glass/epoxy, advances in production techniques and 
developments in the textile industry have resulted in fiber materials with much 
improved strengths and stiffness. Two of the early break-throughs were 
graphite and boron fibers which have typical strengths and stiffness much 
higher than convention ductile isotropic metals. These fibers allow for new 
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design limits but also give rise to the requirement of new design technique 
utilizing a material of an unconventional material redundancy with anisotropy. 

Employing the fiber/matrix as unidirectional lamina to form multidirectional 
laminate gave the designer the tools to make structural designs subjected to 
multiple loading conditions a manageable task. This idea shifted the attention 
of the designer from the microscopic level to the macroscopic level. 

Utilizing the high strength and stiffness, fibers application for structural 
enhancement were incorporated in many of the military designs in the early 
70’s such as the F-lll afl fuselage structure. Testing of these designs were 
accomplished by component. This method is time consuming and expensive but 
must be done. To gain any type of reliability assurance many tests are 
required. It soon becomes apparent that testing costs could out pace production 
costs. 

C. FUTURE DESIGN REQUIREMENTS 

As weaponry becomes more and more sophisticated the demands on 
composite materials for a wide variety of uses are ever increasing. New designs 
often are only limited by the structural limitations of the materials in use. 
Such restrictions might include the thickness of a wing or stabilizer, the sweep 
of the wing, or the "g" loading. Lately new design specifications have included 
the use of composites for the radiation energy absorbing properties. Though the 
Navy’s interest in composites spans a multitude of uses, one special interest is 
in the area of structural enhancement, such as rocket motor cases, pressure 
vessels for submarine flasks and jet aircraft pilot ejection seat. 
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In the last few decades composites have advanced dramatically, and are 
beginning to provide improved alternatives to more conventional methods of 
design normally utilizing isotropic alloys. Fiber reinforced plastics have 
attracted a large amount of attention. Improved techniques in the textile 
industry have lead to more uniform fibers from bases of precursors for high 
performance graphite’s with strengths and stiffness exceeding high yield steels. 

As utilization goes up the increased use of new fibers in composites and the 
cost of system failure testing demand better methods for reliability evaluation. 
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II. BACKGROUND 



A. COMPOSITE LOAD SHARING 

A brittle material generally has a higher strength to density ratio but a 
lack of ductility limits the structural uses. Ductility resists crack propagation 
from a region of damage. Metals have a high ductility and do provide practical 
strength and stiffness properties for many structural applications. Fibers do not 
use ductility to resist failure [Ref. 2:pp. 1-19], because ductility is not required 
the fiber can be synthesized from brittle materials. 

The maximum theoretical elastic strain for carbon hydrogen bonds is 
approximately e mai = 0.1, which is one order of magnitude greater than 
present day strong elastic reinforcing fibers. 

Presently fibers can not obtain strengths any where near the theoretical 
ideal strength of the material but are an order of magnitude greater than 
conventional bulk materials. Fiber strength can also be attributed to the small 
fiber diameter, which limits the flaw sizes. 

Composites are geometrically advantageous to homogeneous material such as 
an all metal or all ceramic material, for its ability to utilize the high strengths 
and stiffness of brittle materials while maintaining strength redundancy through 
the matrix. 

The strength interaction of the fiber and matrix is the key to reliability 
effectiveness of a composite. Consider several fibers in parallel without matrix. 
The load is equally shared by each fiber. If one fiber should fail the failed fiber 
is completely ineffective and no longer supports any load. The load of the 
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remaining fibers increases equally. With a matrix present (Figure 1) a crucial 
phenomenon occurs. At the point of fracture the more ductile matrix absorbs 
the load through shear stresses along adjacent fibers and back to the broken 
fiber itself. Near the fracture site high shear stress develope in the matrix 
along the broken fiber but quickly 
dissipate as the normal stress 
increases. On the two surrounding 
fibers a noticeable stress increase is 
observed adjacent to the fracture. 

The short distance from the broken 
end of the fiber to the point that 
fiber is carrying a full stress is called 
the ineffective length [Ref. 3] 

The ability of the composite to 
adjust to multiple fiber failures is it’s 
strength redundancy, and is 
dependent on the characteristics of 
both the matrix and fiber. There are several fracture mechanisms all of which 
are a function of the bond strength (shear strength of interface), the matrix 
shear strength and the distribution of the fiber strength. The most critical is 
the fiber strength distribution. The distribution of the fiber determines the 
composite reliability. 
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Figure 1 . Local load sharing in a 
Fiber/Matrix composite. 
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B. STATISTICS AND RELIABILITY 

The complexities of composites make the task of analysis virtually 
impossible. New fibers are being developed every year and attempting to 
generate thorough material testing studies as has been done with conventional 
materials is not realistic. But composite reliability can be estimated through 
probability studies. Probability of the composite must be based or inferred from 
accurate statistics of the fiber distribution. 

To determine the reliability of the composite, the strength distribution of 
the fiber must be known. The distribution can be quantified by strength 
statistics (a Empirical Probability Distribution Function (EPDF), or Empirical 
Cumulative Distribution Function (ECDF)). For computational efficiency the 
empirical data can then be modeled with an analytical function and thereby 
characterize the distribution by parameters of the function. To define the 
parameters the distribution function must be identified by model type (i.e., 
Normal, Weibull, etc.). Using the model Probability plot; the closeness of fit of 
the data can be judge and therefore the appropriateness of the model. Selecting 
the model allows for the calculation of the model parameters. 

The relationship between the fiber statistics and the composite probability is 
represented in a linearized Probability plot (Figure 2). 

Reliability is mostly concerned with the lower tail of the distribution, this 
also being the most difficult to define. Because of the redundancy provided by 
the matrix induced load sharing the composite is considerably more reliable 
than it’s constituent fibers in the lower tail. As the fiber statistics fluctuate in 
the extreme lower tail the composite reliability fluctuates with it. This can be 
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Probability Plot 
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Figure 2. Comparison of Fiber and Composite in a 
Probability Plot. 

shown in a representation of a probability plot in Figure 3. There is a critical 
point where the fiber plot stabilizes, the linear region in this plot determines 
the accurate statistics of the lower bound reliability (worse case) of the 
composite. 

Fibers have a characteristic which helps to identify the analytical 
distribution model that will closely define the empirical statistic. Traditionally, 
the weakest link model [Ref. 4] depicted the fiber as being made of many small 
segments linked together much like a chain. The segments each have an 
intrinsic strength, none being exactly equal. When the fiber is stressed the 
weakest will fail and the total fiber fails. The Weibull distribution model 
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defines this weakest link behavior and is commonly used in modeling fiber 
strength. 

1. Weibull Density Distribution Function 

The Two parameter Weibull cumulative distribution function (CDF) is 
defined by the following equation. 

F(x) = 1 - exp {-(x/B) a ), x, B, a >0 (2.1) 

where: x = random variable (stress, load) 

B = scale 
a = shape 

The parameter alpha defines the shape (spread and skew) of the 

distribution. Alpha less than 3.5 the 
plot is skewed positively, greater than 
3.5 it is skewed negatively and 
neutral at 3.5. The alpha parameter 
determines the strength scatter of the 
fiber, an increasing alpha defines less 
scatter and higher strength 
uniformity. The beta parameter 
defines the central tendency of the 
distribution and is relatable to the 
mean. Alpha for a fiber type is 
constant for the weak link physics, 
whereas beta changes with the fiber length. An increase in length causes a 
decrease in beta (a weaker fiber). 
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The composite reliability can be inferred from the distribution of a fiber 
strength. To gain high reliability with high confidence levels a proportionate 
amount of statistics must be obtained (i.e. for a reliability level of 10 3 requires 
10 4 statistics). 

Presently fiber statistics are gathered by single filament testing (Figure 
4). First a single fiber specimen is prepared, since a fiber has such a small 
diameter, (AS-4 diameter ~ 10 micron) it is very fragile. The fiber is placed in 
a cardboard frame. Once the specimen is loaded in the test apparatus the 
cardboard frame is severed by a 
heated wire followed by load 
application and interpretation of the 
load deformation data. This is a 
time consuming effort which 
produces only one failure strength 
statistic. To gain the number of 
statistics necessary would take 
thousands of tests, as an example man-safe reliability is 10 s therefore to 
achieve the necessary statistics would require 10 7 tests. 

A new approach being explored is to test many fibers in parallel (a 
bundle) and recover the fiber characteristics from the failure test. This would 
provide a large statistical base required and could be accumulated with a much 
higher efficiency. With a validated method the analysis of fiber reliability could 
keep pace with the development of new fibers with reduced effort. 




specimen. 
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2. Size Effect. 



The extreme lower tail is difficult to determine accurately for the fiber. 
For a Weibull failure physics and the resulting strength distribution the lower 
tail can be estimated by assuming weakest link [Ref. 4], The parameter beta 
is a function of the fiber length and is defined by equation 2.2 

fi 2 = (2.2) 



where 13, = beta for 1, 

13 2 = beta for 1 2 
a = a, = a 2 

A detailed development can be seen in Appendix A. Using this 
relationship parameters for beta can be retrieved analytically where practical 
testing is impossible and the extreme lower tail is estimated. 

This thesis will attempt to measure the lower tail via bundle testing, as 
opposed to the aforementioned estimation, (fiber statistics via single filament 
testing). This investigation explores the methodology needed in testing a bundle 
of fibers to failure. By recording the load and displacement it is expected that 
the underlying distribution of the fibers in that bundle may be extracted. If, 
for example, a bundle has 3k of fibers, the distribution is based on 3k of 
statistics. As the statistical base approaches the complement of the reliability, 
the accuracy of the estimation increases. 
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III. TECHNICAL ISSUES AND RESOLUTIONS 



A bundle of fibers tested to failure presents several experimental obstacles 
which effects the analysis. There are three experimental consequences to 
testing a bundle in parallel that are not experienced in testing a single fiber in 
series. These phenomenon which may disturb the failure distribution, are 1) a 
protective/adhesive enhancement coating called sizing which is placed on the 
filament spools during production, 2) friction which occurs during testing after a 
portion of the bundle has failed and 3) slack which results from specimen 
fabrication. Each of these is discussed in the following paragraphs. 

A. SIZING 

Fibers are packaged on spools of bundles, in order to prevent the fibers 
from entangling a small amount of a chemical liquid called sizing is coated on 
the bundle. The effect of this sizing on bundle testing is unknown and must be 
address at sometime. The testing in this experiment was done leaving sizing on 
the bundle. A method to remove the sizing chemical without damage to the 
fibers needs to be studied but is not addressed here. 

B. FRICTION 

During the bundle test prior to fiber failure there is negligible friction, this 
is due to the equal strain between fibers. As the fibers fail the strain of the 
failed filaments return to zero leaving the loose ends entangled in the 
remaining fibers, this places added strain on those fibers. The resulting effect 
shows up as reduced failure strengths and clustering of the these failures. In 
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the relatively short gauge lengths friction is expected to be small and very 
pronounced in the long gauge lengths. The effects of lubricating the bundle 
with oil during testing will be observed to see if the effect is reduced. 

C. SLACK 

During specimen preparation and placement in the test apparatus the 
alignment of the fibers becomes misaligned with unequal length. The effective 
result is an initial nonlinear load/displacement which also effect the failure 
distribution. Another concern is the deformation in the actual load train of the 
test apparatus. This is referred to as system compliance. Compliance was 
found to be a function of the load and can be removed from the resulting data. 
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IV. EXPERIMENTATION 



The experimental setup consisted of a Material Tester (INSTRON 4206) 
and an integrated data acquisition system consisting of an Instron Control 
Console, IBM PC/AT and software, a detail description and a list of procedures 
is in Appendix B. The object of the experiment was to test a current composite 
fiber type (used in many Navy aircraft) with a range of several gauge lengths 
and attempt to extract the underlying distribution of the fiber. The sample 
preparation is described in Appendix B. Each bundle sample was tested, the 
data stored, and quick analysis conducted to ensure no improper testing method 
were corrupting the data. 

A. BUNDLE TESTING 

The fiber material tested was made from a Hercules Magnamite high 
strength graphite, type AS-4 spool 145 in 3k bundles (i.e. there are 3,000 
filaments in a bundle). Several different gauge lengths were tested in order to 
analyze the effects of friction, slack, compliance on the fiber size. Size affect 
(Appendix A) should be able to correlate the location parameter beta between 
fiber lengths if weakest link physics is applicable. In all, nineteen bundle 
samples where tested. The test procedure was conducted as specified in 
Appendix B. The complete test log is listed in Table I. 

The bundle samples were handled with extreme care so as not to damage 
the any of the fibers. Moving the bundles from assembly bench to the storage 
bench and then to the testing machine presented the greatest hazard. The 
difficulty came when picking up the sample. This was done manually and 
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Table I. TEST LOG. 






Sample # 


G.L. 

(cm) 


Oil/ 

Dry 


Comments 


090901 


5.0 


dry 


Good 


090902 


0.5 


dry 


No failure-slip 


090903 


0.5 


dry 


No failure-slip 


090904 


0.5 


dry 


No failure-slip 


090905 


0.5 


dry 


No failure-slip 


090906 


0.5 


dry 


No failure-slip 


090907 


5.0 


dry 


Slip 


090908 


5.0 


dry 


Good 


090909 


5.0 


oil 


Good 


090910 


2.5 


dry 


Possible damage 


100901 


2.5 


dry 


Good 


100902 


2.5 


oil 


Slip 


100903 


2.5 


oil 


Slip 


100904 


50.0 


dry 


Good - (friction) 


100905 


50.0 


oil 


Good - (friction 
reduced) 


100906 


50.0 


oil 


Possible damage 


100907 


25.0 


oil 


Good 


100908 


25.0 


dry 


Good 



bending the bundles slightly was unavoidable. The long bundles (500 and 250 
mm) were less vulnerable to damage because of the relatively large aspect ratio 
but these were difficult to steady during transportation and tended to bounce 
around. 

Placing the sample in the grips also presented operational difficulties when 
tightening the upper grip with a socket wrench. An attempt was made to 
manually counter torque the grip, but this tended to be jerky at best. Also the 
grip pressure was difficult to control between samples. The alignment of the 
bundle in the grips was another concern. This was done by sight. If not done 
properly the bundle failure mechanism might be influenced. Manually reducing 
the experimental slack (not to be confused with bundle slack) was very touchy. 
When toggling the grips at the control console the load read out had to be 
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acutely monitored for the bundle loading must be kept to a minimum prior to 
testing. 

During testing a real time graphical display of the load vs displacement 
was used for monitoring the results, with this experimental irregularities in the 
several of the samples was detected. The bundles in the very short gauge 
lengths showed irregular load curves. Keying into this after several samples 
helped to identify a problem with the sample preparation. 

Friction was expected to affect the long gauge length (250 and 500 mm) 
and not affect the shorter gauge lengths (25 and 50 mm). To test the affect 
each gauge length tested a pair was treated with oil and a pair was left 
untreated (dry). Tests were run on samples of 500 mm, 250 mm, 50 mm, 
25 mm and 5 mm. 

B. DATA REDUCTION 

After tests were complete the data files were decoded into ASCII and 
stored for processing. The data output from the Instron software is a three 
column file. The first column is the record number, the second is the 
displacement (inches) and the third is the load (lbm). This is then processed 
through a program to convert the data to mm and kg and remove compliance. 
There are two data files output, one file with compliance removed and the other 
converted raw data (compliance not removed) [Ref. 5]. The output from the 
conversion program is then processed through an analysis program listed in 
Appendix C. This outputs the ECDF, the bundle Modulus E and the Slack 
region ECDF. The ECDF is then processed through a Weibull Maximum 
Likelihood Estimator program to define the shape parameter alpha and the 
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location parameter beta. ECDF is also processed through a linearizing program 
which outputs the Weibull probability plot data file. 
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V. RESULTS 



A. TEST OBSERVATIONS 

During testing several of the specimens were invalid due to experimental 
difficulties. Some of the samples had maximum strains noticeably beyond 
expected limits. One bundle was observed to have a much lower modulus than 
what was expected from equation C.5. Dissection of the samples with suspect 
loading curves revealed a lack of adhesive wetting resulting in a large 
percentage of the fibers slipping from the tab. The problem appeared most 
frequently in the small gauge length (GL) Samples. It is suspected that the 
adhesive which couples the bundle and tab was not wetting the internal fibers 
of the bundle, with a resulting loss of ability to carry the higher loads of the 
small GL. This phenomenon manifested itself two ways. If there was partial 
slipping after initiation of the failure region, a strain higher than the actual 
maximum strain was observed. If total slipping occurred in a large group of 
fibers the modulus decreased. After closer inspection of all the samples it was 
seen that many of the specimens had a small number fibers slipping in the tabs 
but the percentage was small enough not to effect the results. 

To avoid this problem the remaining tests were conducted in samples at or 
beyond 25 mm. The problem was less acute in the long GL because of the 
lower loads these long GL’s experience as a consequence of size effect. 

Two types of test were carried out, with oil and dry. This data was then 
reduced by categories of compliance removed and not removed. The compliance 
was estimated by testing samples of zero GL [Ref. 5]. The resulting curves 
were fitted and reduced to an analytical equation, displacement as a function of 
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load. This equation can then be used to subtract the compliance from the 
empirical data. The zero GL samples used the same coupling adhesive and 
possibly suffered from partial slipping. This causes the compliance curve 
coefficients to be over estimated and as a result displacements removed for 
compliance were too large. When the compliance was removed from the 
experimental data displacements in the lower range showed negative values 
indicating over estimations. 

Since the determination of compliance is slightly in error it is not used in 
the results discussion, though data reduction calculations were completed to 
indicate possible trends and affects. Table II shows the data reduction with 
removal, Table III is without compliance removal and Table IV, for comparison, 
is the result from single fiber tests [Ref. 6]. 

Table II. EXPERIMENTAL DATA RESULTS, COMPLIANCE REMOVED. 



Compliance Removed 



Test # 
dry/oil 


Gage 

length 

(mm) 


Modulus 

(kg/mm) 

(gm/e) 


alpha 


Beta 

(mm/mm) 


100901 

dry 


25 


109.76 


915 


3.0 


0.0157 


090901 

dry 


50 


64.09 


1068 


3.214 


0.0135 


100908 

dry 


250 


11.0 


917 


3.18 


0.00972 


100904 

dry 


500 


5.433 


905 


3.035 


0.0076 


090909 

oil 


50 


54.7 


912 


4.465 


0.0145 


100907 

oil 


250 


11.18 


932 


3.964 


0.0101 


100905 

oil 


500 


5.755 


959 


3.66 


0.0085 



18 



B. DATA COMPARISON 



Generally the results indicate only small deviations in the alpha parameter 
(Table III) for all GL. Between a oiled sample of the same GL only the 500 
mm GL had any noticeable change in the shape parameter, which had been 
expected. The location parameter was effected most by the oil treatment, beta’s 
increased by 10% in both the 50 and 500 mm GL. In 250 mm GL the alpha 
and the beta parameters seem to be relatively stable, meaning the oil treatment 
had little effect. 

It was expected that friction would affect the longer GL’s and have little 
or no affect on the shorter samples. The 500 mm test did show a considerable 
difference when oil was applied (Figure 5). 

The friction was expected to become influential only as the percentage of 
failed fibers increased, but the dry bundle shows a decreased load early, which 
can be seen in the ECDF plot (Figure 6). In fact the oil and dry sample data 
merge as the percentage of failures increases. The reduced friction made a 
difference in the resulting Weibull parameters, alpha increased from 3.085 to 
3.8. From single fiber testing alpha (Table IV) for AS-4 is approximately 4.11. 
Even with oil the bundle is affected, as seen by the sudden drops in the bundle 
load (Figure 5) or the sharp increase in the percentage of failures (Figure 6) at 
discrete points. 

Without friction the failure distributions are expected to be smooth. 
Decreasing the GL it is anticipated that the effects of friction diminishes. Tests 
at 250 mm the loading curve is seen to smooth out (Figure 7), providing a 
slightly more continuous plot. Friction had less effect on this GL though a 
slight decrease is noticed when the bundle is dry. 
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Table III. DATA REDUCTION, COMPLIANCE NOT REMOVED. 



Compliance Not Removed 



Test # 
dry/oil 


Gage 

length 

(mm) 


Modulus 

(kg/mm) 

(gm/e) 


090901 

dry 


50 


45.6 


761 


100908 

drv 


250 


9.85 


820 


100904 

dry 


500 


5.2 


867 


090909 

oil 


50 


43.5 


725 


100907 

oil 


250 


10.6 


833 


100905 

oil 


500 


5.5 


917 



Alpha 


Beta 

(mm/mm) 


Beta 
sized to 
50 mm 


4.75 


0.01615 


0.01615 


4.07 


0.0104 


0.0154 


3.085 


0.00781 


0.0137 


4.89 


0.0179 


0.0179 


3.95 


0.01041 


0.0154 


3.8 


0.00868 


0.0152 



SUMMARY 


a 


Std 


B 


Std. 


E 


Std 


OIL 


4.21 


0.59 


0.0162 


0.0015 


825 


96 


DRY 


3.96 


0.84 


0.0151 


0.00125 


816 


53 


Merged 


4.09 


0.66 


0.0157 


0.00137 


821 


75 



Table IV. RESULTS FROM SINGLE FIBER TEST AT 50 mm GL. 

Fiber Load Weibull Parameters 



Spool # 


Alpha 


Beta 


008 


4.28 


0.0175 


019 


3.94 


0.0172 


average 


4.11 


0.0173 



Modulus = 900 gm/e 
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Figure 5. 500 mm gage length load curve comparison. 



21 



0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 

DISPLACEMENT (mm) 



1 




o o o’ o o 



O) J 



Figure 6. ECDF 500 mm comparison. 
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The ECDF plot (Figure 8) indicates friction had more influence than was 
revealed in the load curve. The upper tails again merge, indicating either the 
oil no longer decreases the friction or friction no longer affects the remaining 
fibers in this region of the ECDF. A slight change can be seen between the 
two samples in the lower and mid region but this could be variations between 
bundles and not whether the specimen was treated with oil, but it should not 
be dismissed because it is consistent with the other samples. 

A 50 mm GL test (Figure 9) showed results which had been anticipated in 
the longer lengths. Little change is observed as few fibers fail but the affects of 
friction become prominent as the percentage increases. 

The ECDF (Figure 10) indicates the effect of oil more clearly. The entire 
range of failure is shifted not only at the upper tail but the lower tail as well. 

Determining how well the underlying Weibull distribution is extracted from 
a bundle test and what influence the oil treatment has can be observed in 
Weibull plots. For the 500 mm sample (Figure 11 & 12) friction has definitely 
affected the distribution, but it is also observed that oil helped to reduce some 
of the deviation especially in the lower tail. 

As the gauge length is reduced it becomes apparent from the Weibull plots 
(Figures 13,14,15 & 16) for the 250 mm GL and the 50 mm GL that friction is 
not necessarily negligible, for the distribution is affected. This is clearly 
demonstrated in the 50 mm GL treated with oil (Figure 15) which shows the 
characteristic Weibull distribution distinctly as apposed to the dry sample. 
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Figure 7. 250 mm gage length load curve comparison. 
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Figure 8. 250 mm Loading curve comparison. 
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'figure 9. 50 mm GL load comparison. 
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Figure 10. ECDF for 50 mm comparison. 
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Figure 11. 500 mm GL Weibull plot, oil treated. 
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Figure 12. 500 mm GL Weibull plot, dry. 
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figure 13. 250 mm GL Wiebull plot, oil treated. 



30 



L0G(X/BETA) 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 m 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1' m rr r r nTi T i 1 1 1 1 1 1 1 1 1 1 1 1 



o 

in 

c\i 

I 



OO'Z 00' L OO'O 00' L- OO’Z- 00’£- 00'f- 00'0- 




Figure 14. 250 mm GL Wiebull plot, dry. 
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Figure 15. 50 mm GL Wiebull plot, oil treated. 
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Figure 16. 50 mm GL Wiebull plot, dry. 
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VI. SUMMARY AND RECOMMENDATIONS 



A. SUMMARY 

The affect of friction is evidenced in all gauge lengths and is not always 
limited to the upper tail. Lubrication of the bundle reduces the effect of friction 
and should be used in all bundle failure testing. 

The compliance removal was not used in the data reduction due to the 
error in the compliance curve caused by adhesive coupling failure in the tabs. 
The compliance ended up being over estimated. This a side the trend of 
compliance removal on the data reduction (Table II) indicated a general 

reduction in the parameters (ot,0). Estimating the compliance to be much less 
and extrapolating this thought, it can be assumed the value for the parameters 
presently determined would be slightly less. In a comparison of the single 
fibers parameters with the bundle parameters (Table III), the bundle 
parameters were barely higher. Given the affect of compliance it can be 

concluded that once removed the bundle parameters will be indistinguishable 
from the single fiber testing and in fact given the quantity of the statistical 
base much more reliable. 

It was shown that friction effected all the gauge lengths but the longer 

gauge length (500 mm) experienced the most disturbance. The 50 mm gauge 

lengths when treated with oil produced a very smooth distribution in the 
Weibull plot. The dry 50 mm sample produced a poor Weibull plot. The 50 
mm gauge length produce good results through out the entire range. The 500 
mm gauge length produced good results in the Weibull plot at the lower tail 
(Figure 7). Since the lower tail of the large gauge length are statistics of much 
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higher reliability than the 50 mm the linear statistics of the 500 mm can be 
pieced into the 50 mm statistical plot there by extending the range of the 
variable (strain). The upper tail of the 50 mm gauge length is not completely 
linear in the Weibull but this can be removed and statistics from a 5 mm test 
can replace it. This piece meal work on the Weibull plot can be done because 
of size effect. The slope alpha does not change but beta decreases in the 
Weibull as the gauge length increases. With this procedure friction and other 
contaminates can be removed and the underlying parameters accurately 
predicted. 

Through simulation a normal or uniform distribution of slack was 
determined to have very little effect on the distribution if the dO displacement 
is subtracted. 

B. RECOMMENDATIONS 

To ensure better results several areas of experimentation can be improved 
on. 

- The handling of the specimens should be reduced, a rigid carrying 
mechanism which can transport the bundle from assembly to storage 
to testing without any human handling or disturbing the fibers is 
necessary. 

- The manual mechanical grips should be replaced with a pneumatic device. 

- Oil treatment should be used on all samples. 

- Improve the tab coupling. 

The slack distribution requires more investigation than was addressed here. 
Simulation was a key element which helped to determine it’s effect on the 
distribution. An analysis algorithm is provided in this thesis which will help to 
determine the slack ECDF. To further study the question of slack several 
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bundle samples (large gauge length) dedicated solely for the determination of 
slack should be tested. A much larger data base is required to use the methods 
discussed in Appendix C. By decreasing the cross-head rate and increasing the 
acquisition rate a large data file can be attained. Once the distribution is fitted 
to a model, simulation can determine what affect it has on the bundle strength 
distribution. 
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vn. CONCLUSIONS 



The results achieved in this thesis conclude convincingly that the 
distribution of the fiber can be retrieved through bundle failure tests. The 
parameters extracted from the testing compared favorably with the single fiber 
testing. As a result, enormous characterization efficiency is realized; in fact 
with this technique, statistics to 10 6 are realized. 
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APPENDIX A 



SIZE EFFECT 

The determination of the scale parameter beta of a gage length, which is too 
short to preclude practical experimental implementation, can be accomplished by 
size effect of the weibull distribution function. 

Starting with the Weibull CDF for a fiber of defined gage length 1,: 



F,(x) = 1 - exp{-(x/B,) Ql } 


(A.l) 


R,(x) = 1 - F,(x) 


(A. 2) 



For a second fiber of gage length 1 2 with, 

1, < 1 2 and 171, = m ; m = integer > 0 
and assuming the reliability R, is constant. 

R, = (R 1 )(R 2 XR 3 )-.(R m ) = R, m (A. 3) 

substituting A. 2 for R, 

R, = expKx/fl.r'P (A.4) 
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and the total reliability 



R t = exp(-(x/13 m ) a ”J (A.5) 

equating equation A.4 and A.5 

expKx/BJH = expMx/B/ 1 }" 1 (A.6) 

taking the natural log A.6 reduces to 

(x/B m r = mlx/B,)" 1 (A. 7) 

again taking the natural log of A.7 

a„ln(x/B m ) - a,ln{m(x/B,)} = 0 (A.8) 

since = a, = a is not = 0 

B,/B 2 = m 1/a (A. 9) 
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APPENDIX B 



EXPERIMENTAL PROCEDURES 

A. EXPERIMENTAL SET UP 

The Test apparatus consisted of an INSTRON Model 4206 material tester, 
with a 50.0 kg load cell, connected to a 4200 Series Expanded Control Console. 
The Console was connected to a IBM PC/AT through com port 1 via a IEEE 
connection and converter. The data acquisition and material tester control was 
commanded by Instron Series IX Automated Material Testing System series 
4.01 C software . 

The controlling software has a 5 page menu which requires setting prior to 
testing. Many of these settings are for Instron data reduction processing which 
are not used in this thesis but required values to operate. The settings crucial 
to testing and proper acquisition are cross-head speed, data acquisition rate and 
gage length. Safety features include a maximum load and displacement setting 
to insure the load cell and testing model are not damaged. 

B. INSTRON CALIBRATION AND TESTING 

The Instron 4206 and control unit required a warm up time of one hour 
prior to operation to allow for stabilization. The load cell has two methods for 
calibration, an electronic calibration and a mechanical calibration. A mechanical 
calibration was performed for all tests. 

1. Mechanical Calibration 

Any inputs refer to the 4200 Series Expanded Control console. 

- Main power on/ofi/on - wait for diagnostic 
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- Press Load Balance / Enter 

- Hand prescribe a 5.0 kg weight 

- Press Load Cal 

- Enter the weight / Enter 

- Remove the weight 

- Enter Load Balance / Enter 

- Re-hang the weight to verify correct calibration 

- Repeat if required 

2. Loading the Sample 

The procedures for sample placement into the test apparatus are as 
follows. 

- The Instron specimen grips are separated enough to 
allow the sample to hang from the top grip without 
touching bottom grip. The tab of the specimen is 
firmly clasped with tweezers and gently lifted off 
resting bench allowing to hang free. 

- The samples top tab is placed in the top grip, 
careful to align the bundle with the center, which 
has been measured and marked. After alignment the 
grip is tightened making sure not to twist or damage 
the specimen. 

- A load balance is performed at the control panel to 
remove the weight of the specimen. 

- The cross-head is then toggle down so as to place 
the bottom tab in the bottom grip then tighten. 

Note: after tightening the grips the specimen develops 
noticeable slack and the load cell measures a negative 
load. The specimen fibers are slightly compressed and 
need to be straightened. 

- The cross-head is toggle up until a load of 0.1 kg 
is indicated in the load readout window on the 
control panel. The load is then brought down to 
zero and the displacement reference zeroized. The 
procedure is repeated once to ensure the 
displacement returns to at or near zero. The 
specimen is ready for testing. 

3. Testing 

The test is initiated at the IBM PC/AT through the Instron software. 
After the specimen test parameters are set and loaded, prior to test initiation. 
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the IEEE port is enabled. The software then signals the system is ready to 
begin testing. 

As the cross-head displaces, load and displacement data are sent to the 
PC/AT at regular intervals of displacement, as prescribed by the initial software 
settings. 

A real time graphical output is available to monitor the testing. This 
allowed for a quick identification of bad tests. After the test was complete, the 
data file is converted from Instron System code to ASCII for latter data 
analysis. 

C. SAMPLE PREPARATION 

The samples were made from a Hercules Magnamite high strength graphite, 
type AS-4, spool 145. The bundle has 3000 fibers with a denier of .0057446 
gm/in.. 

1. Procedures 

A length of bundle is clamped at one end on a 3 meter aluminum track 
and weighted on the other. Half Copper tabs approximately 2.54 cm (width) x 
2.54 cm (length) are pre-positioned on the track at a specified gage length. 
Slots built into the track, center the bundle. An adhesive is applied at the tab, 
the bundle and the remaining half of the copper tab are securely clamped on to 
complete the tab. The bundle is severed between samples on the track and the 
adhesive is allowed to cure. 

Once dry the samples are manually removed from the track and placed 
on the sample bench. 
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D. CALCULATION OF CROSS-HEAD SPEED 



The cross-head calculation is a function of the number of data points 
required, data acquisition rate, load cell and displacement gage tolerance. 

The time to bundle failure was selected to remain constant. This will 
provide a constant strain rate on the different gage lengths that were tested. 

The relationship between the cross-head, acquisition rate and recorded range 
is given in equation B.l. 

XH = (RRXAR)/DPTS (B.l) 

where: XH = cross-head speed 

RR = recorded range 
AR = acquisition rate 
DPTS = number of data points recorded 

The sampling interval is 

deltaD = RR / DPTS = XH / AR (B.2) 

In order not to sample beyond the abilities of the acquisition system the 
following rule must be followed 

deltaD > Recording Tolerance (B.3) 

The Instron output tolerance for displacement was 0.00254 mm, therefore 
deltaD > 0.00254 mm. Using this as a guideline and the constant testing time 
chosen of 5 min.. Table IB lists the Instron test setting calculations. 
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Table IB. TEST SETTINGS. 



Gage 

length 

(mm) 


Max 

displ. 

(mm) 


Time 

(min) 


Cross- 

head 

(mm/s) 


Acq. 

rate 

(pts/s) 


500 


25 


5 


8.33xl0 2 


6.6 


250 


12.5 


5 


4.17xl0 2 


3.3 


50 


2.5 


5 


8.33xl0 3 


3.3 


25 


1.25 


5 


4.17xl0' 3 


1.67 



44 



APPENDIX C 



DATA ANALYSIS 



The analysis of the bundle failure load vs displacement curve can be divided 
into three distinct regions (see Figure 1C): 

- Slack region. 

- Linear region. 

- Non-linear region. 

Slack is created during specimen production and placement in the Material 
testing machine. As the bundle is measured, cut and tabs are placed on each 



LOAD CURVE 



REGIONS 




Figure 1C. Regions of the load curve. 

end of the sample some of the individual fibers become loose and effectively 
longer than other fibers. This produces an uneven loading across the fibers as 
the cross-head of the tensile testing machine strains the bundle. The i lh fiber is 
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not loaded until the displacement of the cross-head is equal to the respective 
slack of the i lh fiber. 

When the cross-head displaces beyond the maximum slack, the loading curve 
becomes linear until the first failure. 

The amount of slack varies according to an undefined distribution. The 
slack distorts the failure region thereby disturbing the underlying sequence of 
the strength distribution. The effect slack has on the underlying shape 
parameter depends on the slack distribution and amount of slack relative to the 
total strain, this is demonstrated in Appendix D. 

A. SLACK DISTRIBUTION 

The slack in a bundle is dependent on the variation in the fiber lengths. 

Lj = 1 + delLj (C.l) 

where: Lj = length of fiber j 

1 = mean gage length 
delLj = slack of fiber j 

Since the delL^ « 1 the length is approximated by 

L, = 1 (C.2) 

The load in the slack region can be defined by Hooke’s law. 

k-1 

P s (d A ) k = I E [(d A - delLp/d + delLj)] (C.3) 

j=i 



k-1 

= I E d/l ; k=l,2,....,n 

j=i 

where: P s ( ) k = bundle load as each fiber k is loaded 

E = fiber modulus (gm/strain) 
d A = apparent displacement = dj + delL, 
dj = fiber displacement 
n = number of fibers 
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Figure 2C. Histogram representation of a PDF plot. 

Defining the loading slack in terms of the distribution function (f s ) of P s at 
a displacement d A as shown in Figure 2C. P s can be write as 

P s (d), = (f s (d A1 )(d A2 -d A1 )n) E U(d-d A1 ){(d-d A] )/l) (C.4) 

where: d = any displacement such that d > d A , 

U(d) = unit step function {=1 d > d A1 ; 

{=0 d < d A1 

Summing the load (eqn. C.4) over the range of d A and defining the bundle 
modulus as a function of the fiber modulus 

E = nE/1 (C.5) 

P s (d) = E I [f s (d A jXd A UI - d A *)] U(d-d A ■) (C.6) 
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differentiating with respect to d for a continuous function 



fd, 

8P/8(d) = E f s (d) d(d) = E F(d_) (C.7) 

Jd, 

Fffip is the CDF of the slack. Assuming the empirical data approximates 
the continuous data. 

P’ = F = F = P’ (C.8) 

The derivative of the empirical data can be numerically solved with two 
methods, a numerical differentiation or a discretization. 

1. Numerical Differentiation 

Even interval Forward and Central difference methods are used due to 
the nature of the slack loading and the recording system. As the fibers load 
the slope is discontinuous as seen in Figure 3C. Several numerical methods 
exist which will produce varying 
accuracy. In trying to recover the 
underlying distribution of the slack 
the sampling rate and number of 
fibers must be considered (i.e., how 
many data points can be recorded 
relative to the number of fibers). 

If the sampling rate is fast 
and the number of total data points i 
recorded in the slack region is 2-3 
times that of the fiber, a 3 point 




• underlyini loaiinfl 
D discrtle data 



Figure 3C. Slack region showing 
discrete data on the underlying loading 
slope. 
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forward difference method may yield accurate results due to the linearity 
between discrete data. If the relative sampling rate is slow the 3 Pt. forward 
difference may skew the distribution. To prevent this a central difference 
method is safe. 

3 pt. Central Difference method. 

P- = (-Ei., + Pi„)/2h (C.9.1) 

5 pt. Central Difference method. 

P ’ = (P , 2 -8P , r 8P U1 -P U2 )/2h (C.9.2) 

where: h = d A ui - sU 

2. Discretization 

The data in the slack region can be fit to a curve and a equation 
estimated. This equation can be differentiated analytically. The underlying 
data itself is not smooth so such a method is not recommended. 

B. LINEAR REGION 

The slope of the linear region is the bundle modulus. To calculate this 
slope a least squares method is very accurate, due in part to the natural 
linearity of the data. Before this method can be employed the region must be 
specifically defined. The point where the non-linear slack region ends is 
identified as dl and the end of the linear region is identified as d2 (the point of 
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the first failure). This is graphically displayed in Figure 4C. Between dl and 
d2 the least squares method calculates modulus E. 

The points dl and d2 must be graphically identified and manually entered 
into the analysis program. To do this a plot of the load curve must be 
generated before the analysis can be conducted. 




1. Least Squares Method 

Due to the noise, the discrete data i will deviate slightly from the 
underlying slope E. From Hooke’s law the expected relationship is 

P, = = ED, (C.10) 

where: D* = d A fdl (C.10.1) 
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Figure 5C. Least Squares Method. 

But due to noise the load P deviates by an amount Pe, this is 
graphically represented in Figure 5. 

P^ = ED, - Pi (C.ll) 

Defining the least squares sum M 

M(E) = I {(ED, - P^o 2 } 2 (C.12) 

where: a = standard deviation 

Differentiating M with respect to E and setting it equal zero defines the 
minimum least squares. 

dM/dE = 0 
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substituting in C.10.1 this reduces to 



E= IP, (tL-dl )/£(&, -dl) 2 (C.13) 

2. Intercept of the Abscissas 

With the Modulus E the intercept of the abscissas can be calculated. 

dO = dl - P(dl)/E (C.14) 

C. FAILURE REGION 

This region is the key to retrieving the strength distribution. Using the 
slope of the bundle the failure statistics can be transferred to a ECDF space. 
From the ECDF the parameters of the distribution can be estimated. 

The ECDF is a plot of the percentage of failures versus a variable. 
Typically this variable is the strength, this would be difficult to translate due to 
the nature of the bundle test. The displacement is constantly incremented, 
where as the load is a bundle load and can vary. Because of this the ECDF is 
plotted versus the displacement. The displacement is proportional to the 

strength as defined by Hooke’s law and once the fiber modulus is identified the 
displacement can be translated to load. 

When a fiber fails, the bundle modulus E changes by one fiber strength. 
From Hooke’s law it is expected that the failures would occur along the 
modulus Ej. Ej being the modulus of the unfailed bundle which ideally is a 
multiple of the single fiber modulus proportional to the number of remaining 
fibers. Unfortunately this is not the case, as the bundle loads, friction and 
slack change the distribution causing the actual loading path to deviate from E, 
and follow along E ( . 
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To show this in more detail consider first the ideal underlying loading (no 
slack, friction or noise) of a bundle under constant strain. At a break in one 
fiber j, the load will drop until intercepting the slope of E,., shown in Figure 6C. 
This new slope emanates radially out from the intercept of the abscissas, (note: 
without slack the intercept should be the origin). If the data acquisition rate 
where high enough this failure could be clearly identified. Now add slack, this 
extends the failure 
displacement of some 
fibers and the intercept 
of the abscissas by E, 
after the initial failure 
to some point between 
the origin and the dO 
point. The intercept 
will shift a relatively 
small amount but wall 
digress toward the 
origin. If friction is 
present the broken fibers will cling to the unbroken fibers and cause additional 
strain. Add noise to the acquisition system and the discrete data fluctuates. 
These factors combine to force the actual recorded data to deviate from the 
underlying bundle loading path. 

A bundle with n number of fibers will have n failures. If the data 
acquisition system can not record this many points in the failure displacement 
range than each statistic can not be recorded. With slack, friction and noise 
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even with a fast enough data acquisition system the failures might not be 
distinguishable. 

The percentage of failures in the bundle test can be defined by the ratio of 
the modulus of no failure to the modulus of failures. The reliability is defined 
by 

& = E/E = ECdJ/tECd^-dO)] (C.15) 

The percent failure is 1 minus reliability, the ECDF is defined by 

F(iJ = 1 - E, (C.16) 

This relationship is graphically depicted in Figure 6C. The ECDF begins at 



Load 

P 



ECDF 
F (%) 




Figure 7C. Transfer from loading curve to ECDF space. 
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d2, so the manual selection of this point must be done with care. 

Assuming all conditions ideal the ECDF would fall out as 

F = j/n (C.17) 

But because conditions are not ideal this does not occur. What is important is 
the shape of the ECDF and whether it is close to the underlying CDF. The 
ability to distinguish each failure is unimportant. Another point to note is the 
lower tail needs to be clearly defined. Due to the relatively low failure rate in 
the lower tail and a constant recording rate this area of the curve has a high 
number of discrete data points per failure, so in essence the testing method 
takes care of itself. 

1. Upper Bound 

The data is recorded on even intervals of displacement not just at the 
failure points, this is graphically described in Figure 6. The discrete data most 
closely resembling the failure point is on the top edge of the loading curve. This 
is referred to as the Upper bound, 
space the Upper bound is on the far 
edge of the plot as seen in Figure 8C. 

To retrieve the underlying CDF from 
the ECDF the Upper bound (Fu) 
must be extracted. 

Attempting to extract the 
Upper bound from a bundle of n 
fibers the most failure points that 
could be optimistically retrieved is n. 

With a large n this is not possible 
due to the acquisition rate, no 



When this is translated into the ECDF 




Figure 8C. Close view of the ECDF. 
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matter, the shape and location is important and does not require a large 
number of points. 

Using the empirical ranking as a reference point 



where m is the number of points desired in the ECDF. F is used to define 
search bands of data where the upper and lower band is defined by 1/n. The 
value of F, which has the largest displacement (£1*,) is assigned to Fu h . When 
no discrete data is in the search band F, Fu k is assigned the value of F at a 
displacement d A k l . 

2. Maximum Likelihood Estimator. 

With the Upper bound extracted this data can be used to estimate the 
parameters of the distribution. The model selected is a Two Parameter Weibull 
Distribution. The estimator chosen is the Maximum Likelihood Estimator 
(MLE) [Ref. #:p. 103] developed by R. A. Fisher. This method finds the 
maximum of the likelihood of the sample in the range of the variable dn and is 
a function of distribution parameters. For the Weibull model the method solves 
for the two parameters alpha and beta. The parameter beta is solved explicitly 
[Ref. 4] in terms of alpha. Alpha can not be solved explicitly and an iteration 
must be used. 



F = k/m ; k = 1,2 ,m 



(C.18) 



B = [ 1/m I du , 2 ] 1/2 sum i=l to m 



(C.19) 



a = { I dn, 2 ln(du,)/(E du, 2 ) - I lnCdu,)}' 1 



(C.20) 
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An Atkinson method of iteration was chosen, initially alpha is assumed. From 
this successive alpha’s are calculated using equation C.20. 



a’ 0 = a 0 - (a 0 - a,) 2 /(a 0 + 2a, + ocj) (C.21) 

where: a„ = initial start value 

a, = f o (a 0 ) 

Ch - f«(a,) 

This equation is iterated for a’,, until a’ 0 = a 0 . 

1. Linearized ECDF 

The data is linearized for a Weibull probability plot to visually check 
the conformity to this model. To linearize the ECDF start with the Weibull 
equation 



F = 1 - exp(-P/fi)“ 

subtracting 1, multi plying by -1 and taking the natural log. 

ln(l-F) = (-P/J3)° 
again taking the natural log 

ln[-ln( 1-F)] = cdn(P/fl) (C.22) 

C. 22 is a linear line with a slope of alpha. 

D. DATA OUTPUT 

As the analysis program calculates certain data file are output, these files 

are: 

• EXPER.OUT, the converted experimental output displacement, strain 
and load. 

- SLACK.OUT, the slack region. 

- EXMS.OUT, the load curve minus slack region. 
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UBFECDF.OUT, the Upper bound ECDF. 
FECDF.OUT, the ECDF. 

SECDF.OUT, slack ECDF. 
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APPENDIX D 



A. SIMULATION 

The simulation of a bundle test has two purposes. First by analyzing in 
depth the mechanisms which govern the failure of the fibers in a bundle, a 
better understanding of how to extract the characteristic properties is attained. 
The properties desired are the distribution parameters. 

Secondly a simulation gives us the ability to know what the underlying 
parameters are and thereby proof the analysis data reduction algorithm. Also 
the variables which influence the shape of the failure distribution can be varied 
in simulation and the degree of induced affect observed. 

The bundle loading curve can be divided into three regions; 

- Slack 

- Linear 

- Failure 

Initially the simulation requires input of the variable estimates, this is 
done through an input file. 

1. Initial Inputs 

The underlying parameters define the strength distribution for a gage 
length. The distribution chosen is the Weibull model. The parameters are 
alpha and beta, the number of fibers per bundle is n. The fiber modulus (E) is 
in gm/mm 2 , from Hooke’s law 

a = E e (D.l) 

where: a = stress 

e = strain 
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Substituting the definition of stress and strain for fiber i 



P/A = E 5/L, (D.2) 

where: P = load 

A = area = constant 
L = length of fiber 
5 = change in length L 

rearranging D.2 and letting E f = EA (gm/e) 

P, = E r 5/L, (D.3) 

2. Slack and Strength Distributions 

The slack is caused by the variation in fiber length L. If L were 
constant no slack would be evident. The slack is expected to have a definable 
distribution, but this is unknown and could vary from bundle to bundle 
therefore several models are simulated. Two models selected where 

- Two parameter Weibull 

- Uniform. 

Although physical consideration favors the slack model to be normal, 
the Weibull model was chosen for it’s analytical convenience. The Weibull is 
very flexible and can characterize a wide range of shapes. For an alpha of 3.5 
the distribution is approximately normal, other than 3.5 the shape of 
distribution is skewed negative or positive. Beta identifies the central tendency 
of the distribution, its relatable to the mean. From equation 2.1 the Weibull 
CDF is defined as 



F w (x;a,fi) = 1 - exp[-(x/B) a 
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rearranging in terms of the variable x,, which can represent slack (delL,) or 
strength P s . 



x, = exp{ [ln(-ln(l-F w ) + aln(B)] / a) (D.4) 

0 < F w < 1 

The uniform model is a simple two variable distribution with a lower 
(u,) and upper (u 2 ) range. Since the slack is initiated at zero the lower range is 
zero. 



F u (delL;u,,u 2 ) = (delL - u,)/(u 2 - u,) (D.5) 

0 < F u < 1 

solving for delL, from D.5 and setting u, = 0 

delL, = u 2 F u (D.6) 

3. Data Generation and Continuity 

A uniform random number generator creates an array of numbers 
between 0 and 1 and stores it until recalled. 

Recalling equation D.3 and substituting in C.l and solving for the 
displacement (5,). 



8, = P, (1 + delL,) / E r (D.7) 

The bundle is put under a constant strain rate, so the discrete 
variable will be displacement. As the bundle is strained the filaments will fail 
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in order of the underlying strength or strain (neglecting slack and friction). For 
a constant filament length the strain is proportional to the displacement and 
the fiber will fail at some underlying displacement. If slack is present in the 
bundle each fiber will not fail at the underlying displacement because the 
displacement is no longer proportional to strain but a function of both slack and 
stain. 

The displacement observed by the Instron material testing machine is 
not 8. The machine records an apparent displacement (d A ) 

d A , = 8, + delL, (D.8) 

This will affect the fiber failure order, every other variable is then 
keyed to the order (r^) of d Aj . 

4. Slack Loading 

As the bundle is strained, the fibers are loaded by a progressive 
summation of each fiber in order of their slack distribution. The loading 
function in the slack region can be defined as a function of the displacement (8.) 
each fiber experiences after slack is released. 

H 

P.,= £E ( S,/L, ; j= 1,2, ...„n (D.9) 

i=l 

j-1 

P §j = 2 E, [(delLo j - delLo,)/(l + delL 0i )] (D.10) 

i=l 

where: delLo , = order slack 
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This relationship can be represented graphically in Figure ID. The last load 
calculated from equation D.10 will occur at delLo n . 




Figure ID. Slack region load simulation. 



4. Failure loading 

The order of each fiber failure is determined by the ordered strength 
distribution (assuming no slack). With slack the order of failure will deviate 
slightly. The failure load is a function of the modulus of the bundle and the 
displacement of each fiber. The bundle failure load (P u ) can be defined from 
equation D.2 as 



P.j = £ E r 8/4 

H 

substituting in C.l and D.8 



P uj = I E [(d AJ - delL,)/(l + delLj 

i=J 



(D.ll) 



(D.12) 
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5. Discrete Data 

The material testing machine measures and transmits data on even 
intervals of displacement. To simulate this the defined loading points 
(P.^elLo) and (P u ,d A ) will be set limits for a linear interpolation. The data is 
interpolated on intervals of deld A which is set by the cross head speed. 
Defining the discrete apparent displacement 

= cL i + deld A (D.13) 

The interpolation algorithm to define the discrete data from 
continuous is the same for both slack and linear region. 

The discrete load P in the slack region is a function of P,, the ordered 
slack (delL 0 ) and the apparent discrete displacement 

E = • P tJ )/(delL 0j ., - delL 0j )]d Ai (D.14) 

The discrete data is defined this way until the first underlying 
displacement d A *. 

The slope of the linear region intercepts the abscissas at a point dO. 
The slope of the linear region is the bundle modulus E (see Figure 2C). 

E = (P„ , - P. „)/(d A , - delLJ (D.15) 

The slope can be solve using E 

dO = delE - P, JE (D.16) 
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LOAD CURVE 




Figure 2C. Linear calculations. 

Subsequent discrete data points can be created by a regressive 
process. At each fiber loss a new slope emanates radial outward from dO. 
Between breaks, to define the discrete data, the slope is defined by the failure 
load P u , the displacement d A and dO 

E = [P u /(d Aj - dO)] (sL, -do) (D.17) 

where: d A , d A , ; j=l,2 n 

6. Noise 

Noise was simulated simply by setting a tolerance level and using the 
random number generator to produce a set of random numbers between -0.5 
and 0.5. This data was added to each discrete load P. 
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The input file to the simulation has units of kg, mm and the 
experimental output from the Instron Software is in lbm,in. A conversion of the 
data is required, this is accomplished through a separate program 

B. SIMULATION INTERPRETATION 

The simulation was written in fortran code, listed in Appendix F. Multiple 
simulations were run varying the parameters slack and noise to see the effect. 

Simulation data was analyzed via the algorithm as represented in 
Appendix C. The output is the bundle modulus E, the shape parameter of 
distribution Alpha, and the location parameter Beta (mm/mm). 

The input parameters were alpha = 4.0, beta = 0.16 kg and a 50.0 mm 
Gage length (GL). The noise was varied from zero to about 0.1% fluctuation of 
the maximum load. The maximum slack input was 3% of the gage length and 
a normal slack distribution. The seed used for the random number generator 
was held constant so the underlying strength distribution was constant. The 
seed for the slack region and noise was not held constant. The results of 
multiple simulation test is in Table IC. A value of 4.134 for the shape 
parameter was retrieved consistently (same seed), as expect no deviation for the 
same data. The increase in slack drove the parameter down. The noise did not 
have any effect on the parameters. Figure 3C show the loading curve with 
slack (minus dO) and no slack. 

No noticeable change in the parameter output nor in the ECDF plot 
(Figure 4C). 
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Table IC. SIMULATION OBSERVATIONS AND 
PROOF TESTING RESULTS. 



Maximum Slack 


alpha 


beta 


o% GL 
(B = o/2) 


(kg) 




0 


4.134 


0.016 


5 


4.010 


0.016 


10 


4.099 

4.095 

4.091 

4.089 

4.088 

4.075 

4.076 


0.016 



The numerical derivative of the slack region is very sensitive to the any 
noise. This can be seen in the plot of the slack ECDF Figure 5C. In order to 
retrieve the slack distribution the recording system must have a low tolerance 
load cell. 
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figure 3C. Simulation load curve, slack; no slack. 
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Figure 4C. The simulation ECDF for slack and no slack. 
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SIMULATION 
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APPENDIX E 



I. DATA ANALYSIS CODE 
A. ANALYSIS PROGRAM 

1. Data Input 

The program requires a formatted data file called EXPER.OUT. The 
data is read in and prompts the user for several variables. These are fisted 
below. 

- What is the gage length? 

- How many fibers in the bundle? 

- Which Numerical differentiation is desired? 

3,5,7 pt. forward diff. or 3,5 pt. central diff. 

- The end of the slack region and the beginning of the 
failure region in displacement (mm) (this requires a graphical 
estimate). 

- How many points in the input data file? 

- The acquisition tolerance of the displacement? (this is for 
any filtering of data points in the event the sampling rate 
is too high and repeated data is evidenced. 

- The modulus is calculated and displayed and then ask if this 
value needs to be change. This is in the event the modulus 
calculation in correct a manual entry can be made. (Be Careful 
not to corrupt the true data). 

2. Data Output 

The data out put is to the screen and to a file. The following is a 
fist of the file outputs. 

* EXPER.OUT, the input file is output with an added column of 
strain 

- SLACK.OUT, the slack region (load vs displacement) 
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- EXMS.OUT, the load curve minus slack region 

- UBFECDF.OUT, the Upper bound ECDF 

- FECDF.OUT, the ECDF 

- SECDF.OUT, the slack region ECDF 

B. MAXIMUM LIKELIHOOD ESTIMATOR 

The MLE program reads one of the ECDF data files, users choice, the first 
line in the input file is a utility line defining the file size, the number of 
filaments and the bundle modulus E> the remaining lines are data. 

1. Data Input 

The input other than the read file is prompted from the screen. The 
following is a list of requested information. 

- How many points for the Expected Ranking (The MLE requires a 
probability ranking) 

- First approximation for alpha 

2. Data Output 

The parameter estimations alpha and beta (mm/mm) are output to the 

screen. 

C. LINEARIZED PLOT 

This program takes the ECDF and does a Weibull linearizing (equation 
C.20). 
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II. FORTRAN CODE 



A. ANALYSIS PROGRAM 
PROGRAM ANALY 

* BUNDLE FAILURE ANALYSIS 



* Outline of program: 

* Read in the data points 



* NLA - Number In Array 

************* ************************ ******** ***************** 



Find the 

-linear region 

-slope of the linear region = EB (E BAR) 

-intercept point of the linear region with 
the abscissas = DO (d zero) 

-Maximum slack or minimum linear displacement - DPI 
-Maximum displacement of the linear region - DP2 






* Input the range of the linear region by manually retrieving 

* the defining points DPI and DP2. 

* 

* OUTPUT: SLACK.OUT - slack region. 

* SECDF.OUT - slack ECDF. 

* EXMS.OUT - experimental data minus slack region. 

* FECDF.OUT - failure region ECDF. 

* UBFECDF.OUT - upper bound of the failure region ECDF. 

* 

INTEGER NIA,NISA 
PARAMETER(NIA=6000NISA=1000) 

REAL D(NIA),P(NIA),PP(NIA),H,LP(NIA),LD(NIA),FS(NIA),F(NIA) 
+,FUB(NIA),DF(N1A),DFUB(NIA),PF(NIA),PFUB(NIA),TEMP(NIA),GL 
INTEGER COUNT, NIDP1,NIDP2,NI(NIA),SLCT 
CHARACTER*!, Q1,Q2 

* 

***************************************************************** 

* Input analysis control parameters. 

* 

***************************************************************** 

PRINT *, What is the gage length? (mm)’ 

READ *, GL 

* 

PRINT *, ’ How many filements?’ 

READ *, N 

PRINT *, ’ Which method do you require for reduction of 
PRINT *, ’ the slack region to the ecdf space?’ 
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* 



PRINT *, ’ 
PRINT *, ’ 
PRINT *, ’ 
PRINT *, ’ 
PRINT *, ’ 
READ *, SLOT 



1) 

2 ) 

3) 

4) 

5) 



3 PT forward diff. method.’ 
5 PT forward diff. method.’ 
7 FT forward diff. method.’ 
3 FT central diff. method.’ 
5 FT central diff. method.’ 



* * * * * * * * * * * * * * * * * * * * * * 4c * * * * * * * * * * 4c 4c 4c * * * * * * * * * * * * * * * * * * * * * * % * * * * * * 

* Read the range of the linear displacement. 

* 



******* **************************************** ****************** 

* 

PRINT *, ’INPUT DPI (the maximum slack)’ 

READ *, DPI 

PRINT *, ’INPUT DP2 (the maximum linear displacement)’ 
READ *, DP2 

******************************************************************** 

* Find the how many points are in the data file. 

* 

******************************************************************** 

PRINT*, ’How many points in the data file?’ 

READ *, COUNT 

******************************************************************** 

* Find out what the tolerance of the machine is. 

* (displacement) 

******************************************************************** 



PRINT *, What is the tolerance of the displacement output’ 
PRINT *, ’data? This is the limit between data displacement’ 
PRINT *, ’Input in mm.’ 

READ *, TOLD 

PRINT *, TIow many points in the empirical rank?’ 

READ * NER 

144 PRINT *, T)o you want the load in Kg (K) or Newtons (N)?’ 

REACH*, 303) Q2 

IF(Q2 .EQ. TC .OR. Q2 .EQ. ’N’) THEN 
PRINT *, ’ ’ 

ELSE 

PRINT *, Try again in UPPER CASE KORN 5 
GOTO 144 
END IF 

* 

******************************************************************** 

* Read in the experimental results. 

* 

******* ******** ***************************************** ************ 



OPEN(UNTT=ll,FILE=’EXPER.OUT,STATUS=’OLD’) 

* 

DO 10 1=1, COUNT 
READU1.110) NI(I),TEMP(I),P(I) 

10 CONTINUE 
REWINDUl) 
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CONVER = 1.0 
IF(Q2 .EQ. ’N’) THEN 
CONVER = 9.8062 
END IF 

DO 44 1=1, COUNT 
P(I) = P(I)*CONVER 

WRITE(1 1,1001) NI(I),TEMP(I),P(I),TEMP(I)/GL 
44 CONTINUE 

************** * ** *************************************** ************ 

* Filter out any points with difference in displacement interval 

* less than the machine tolerance TOLD. 

* 

******* ******* ****************************************** ************ 
J=1 

D(l) = TEMP(l) 

DO 11 1=2, COUNT 
DIFF = TEMP(I)-TEMP(I-1) 

IFCDIFF .LT. TOLD) THEN 
PRINT *, ’SKIPPED POINT ’,1 
ELSE 
J=J+1 

D(J) = TEMP(I) 

P(J) = P(I) 

END IF 

11 CONTINUE 
COUNT = J 

* 

******** 

* FIND COUNT (NI) OF DPI AND DP2 
******** 



1=1 

20 IF(D(I) .LT. DPI) THEN 

1 = 1+1 
GOTO 20 
ELSE 
NIDP1 = I 
DPI = D(I) 

ENDIF 

30 IF(D(I) .LT. DP2) THEN 

1 = 1+1 
GOTO 30 
ELSE 

NIDP2 = 1-1 
DP2 = D(I-1) 

ENDIF 

* 

******************************************************************** 

* OUTPUT THE SLACK REGION. 

* 

******************************************************************** 
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0PEN(UNIT=13,FILE=’SLACK.0UT,STATUS=UNKN0WN) 

DO 500 I = 1,NEDP1 
WRITE(13,110) I,D(I),P(I) 

500 CONTINUE 
CLOSEU3) 

* 

******************************************************************** 

* LEAST SQUARES TO DETERMINE THE SLOPE BETWEEN DPI AND 

DP2. 

* 

******************************************************************** 



TOP = 0.0 
BOT = 0.0 

DO 40 I=(NIDP1+1),NIE>P2 
LP(I) = P(I) - P(NIDPl) 
LEKI) = IXI) - DPI 
TOP = TOP + LP(I)*LD(I) 
BOT = BOT + LEKI)* *2 



40 CONTINUE 

EB = TOP/BOT 
* * * * 


* 


* 


* 


** 


* Find the intercept point do 

* * * * 


* 


* 


* 


** 



DO = DPI - P(NIDP1)/EB 
PRINT *, ’Do = ’,D0 
PRINT *, ’ ’ 

PRINT *, ’The Modulus calculated from the load curve’ 
PRINT *, ’is \EB 
PRINT *, ’ ’ 

PRINT *, ’Would you like to change this (Y/N)?’ 
READ(*,303) Ql 

IF(Q1 .EQ. 'Y’ .OR. Ql .EQ. ’y’) THEN 
PRINT *, ’Input the new modulus’ 

READ *, EB 
ENDIF 

303 FORMATXAl) 



* OUTPUT THE EXPERIMENTAL DATA MENUS SLACK 

* 

******************************************************************** 

OPEN(UNIT=13,FILE=’EXMS.OUT\STATUS=’UNKNOWN) 

J=1 

WRITE(13,110) J, 0.0, 0.0 
DO 550 I = NIDP1, COUNT 
«J=<J+ 1 

WRITE(13,110) J,(D(I)-D0),P(I) 

550 CONTINUE 
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********* ************************************************* ********** 

* NUMERICAL DIFFERENTIATION OF THE SLACK REGION 

* (using a forward and central difference methods). 

* 

******************************************************************** 

♦ 

IF(SLCT .EQ. 1) GOTO 103 
IF(SLCT .EQ. 2) GOTO 105 
IF(SLCT .EQ. 3) GOTO 107 
IF(SLCT .EQ. 4) GOTO 130 
IF(SLCT .EQ. 5) GOTO 150 

* 

********************************** 

* 3 pt forward diff. method. 

********************************** 

* 

103 DO 3 I=1,NIDP1 

* if out of the slack region leave the loop 

IF (D(I) .GT. DPI) GOTO 51 
H = ABS(D(I+1)-D(I)) 

* 3 point forward difference numerical differentiation. 

* Generate Fs(da) CDF 

PP(I) = (-3*P(I)+4*P(I+l)-P(I+2))/(2*H) 

FS(I) = PP(I)/EB 
3 CONTINUE 
GOTO 51 

********************************** 

* 5 pt forward difference method. 

********************************** 

105 DO 5 I=1,NIDP1 

* if out of the slack region leave the loop 

IF (D(I) .GT. DPI) GOTO 51 
H = ABS(D(I+1)-D(I)) 

* 5 point forward difference numerical differentiation. 

* Generate Fs(da) CDF 

PPG) = (-25*P(I)+48*P(I+l)-36*P(I+2)+16*P(I+3)- 
C3*P(I+4))/(12*H) 

FS(I) = PP(I)/EB 

5 CONTINUE 
********************************** 

* 7 pt forward difference method 
********************************** 

107 DO 7 I=1,NIDP1 

* if out of the slack region leave the loop 

IF (D(I) .GT. DPI) GOTO 51 
H = ABS(D(I+1)-D(I)) 

* 7 point forward difference numerical differentiation. 

* Generate Fs(da) CDF 

PPG) = (-147*P(I)+360*P(I+l)-450*P(I+2)+ 
+400*P(I+3)-225*P(I+4)+72*P(I+5)-10*P(I+6))/(60*H) 

FS(I) = PP(I)/EB 
7 CONTINUE 
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********************************** 

* 3 pt central difference method 
********************************** 

* 

130 1=1 

H = ABS(D(I+1)-D(I)) 

PP(I) = (-3*P(I)+4*P(I+l)-P(I+3))/(2*H) 

FS(I) = PP(I)/EB 
DO 33 I=2,NIDP1 

* if out of the slack region leave the loop 

IF (EKI) .GT. DPI) GOTO 51 
H = ABS(D(I+1)-D(I» 

PPG) = (-P(I-1)+P(I+1))/(2*H) 

* Generate Fs(da) CDF 

FS(I) = PP(I)/EB 

33 CONTINUE 
********************************** 

* 5 pt central difference method 
********************************** 

* 

150 1=1 

* Use a 3 pt forward difference method for first couple pts. 

H = ABSOXI+D-DG)) 

PPG) = (-3*P(I)+4*P(I+l)-P(I+3))/(2*H) 

FS(I) = PP(I)/EB 
1=2 

H = ABS(D(I+1)-DG)) 

PPG) = (-3*PG)+4*PG+l)-PG+3))/(2*H) 

FS(I) = PP(I)/EB 
DO 55 I=3,NIDP1 

* if out of the slack region leave the loop 

IF (D(I) .GT. DPI) GOTO 51 
H = ABS(D(I+1)-DG)) 

PPG) = (PG-2)-8*PG-l)+8*PG+l)-PG+2))/(12*H) 

* Generate Fs(da) CDF 

FS(I) = PP(I)/EB 
55 CONTINUE 

51 CONTINUE 

******************************************************************** 

* FAILURE REGION 

* 

******************************************************************** 

J=0 

DO 60 I=NIDP2, COUNT 
J=J+1 

F(J) = 1 • PG)/(EB*(DG)-D0)) 

DF(J) = D(I)-D0 
60 CONTINUE 
NUMF = J 

******************************************************************** 

* EXTRACT THE UPPER BOUND. 
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* 

******************** **************************************** ******** 



DEL = 1.0/NER 

K=NUMF 

I=NER 

FUB(I+1) = F(K) 

DFUB(I+1) = DF(K) 

K=K-1 

1 = 1-1 

FL = 1.0*I/NER - DEL/2 
FU = 1.0*I/NER + DEL/2 
PRINT */FL F(K) FU’ 

PRINT *, FL/ < \F(K)/ < \FU 

90 IF((K .LT. 1) .OR. (I .LT. 0)) GOTO 92 

91 IF((F(K) .LT. FU) .AND. (F(K) .GE. FL)) THEN 

FUBG+l) = F(K) 

DFUBG+l) = DF(K) 

1 = 1-1 

FL = 1.0*I/NER - DEL/2 
FU = 1.0*I/NER + DEL/2 
K-K-l 

IF(K .LT. 1) GOTO 92 
GOTO 91 

ELSEIF(F(K) .LT. FL) THEN 

******** 

* check to see how many it skips 

******** 

J1 = 1 

93 IF(F(K) .LT. (1.0*(I-J1)/NER - DEL/2)) THEN 

J1=J1+1 
GOTO 93 
ENDIF 

DO 94 I1=I,(I-J1+1),-1 
IF(I1+1 .LT. 1) GOTO 92 
FUBGl+l) = 1.0*I1/NER 
DFUB(I1+1) = DF(K) 

94 CONTINUE 

1=11 

FUB(I+1)=F(K) 

DFUB(I+ 1 )=DF(K) 

1 = 1-1 

FL = 1.0*I/NER - DEL/2 
FU = 1.0*I/NER + DEIV2 
K=K-1 
GOTO 90 
ELSE 
K=K-1 

IF(K .LT. 1) GOTO 92 
GOTO 90 
ENDIF 

92 CONTINUE 

PRINT *, ’OUT OF LOOP WITH K AND I ’ 
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PRINT *, ’K= ’,K,’ 1= M 
NUMEF = 1+2 

*** * * *********************************** *********** + *++ 4 :****** 

* WRITE THE DATA 

* 

** ** *** * ** * * ***** ** * * ** ** * *** *** *********** ****** 3t= * ** *4= ****** * 

OPEN (UNIT=1 l,FILE=TJBFECDF.OUT’,STATUS=TJNKNOWN’) 
MSUM = NER+1 -NUMEF 
WRITE(U.IOIO) MSUM,N,EB 

* 

DO 250 I=NUMEF,NER 
WRITE( 11,110) I,DFUB(I),FUB(I) 

250 CONTINUE 
CLOSE(ll) 

OPEN(UNIT=ll,FILE=’FECDF.OUT’,STATUS=TJNKNOWN’) 
WRITE! 11 ,1010) NUMF,N,EB 
DO 70 I=1,NUMF 
WRITE(11,110) I,DF(I),F(I) 

70 CONTINUE 
CLOSE! 11) 

OPEN(UNIT=12,FILE=’SECDF.OUT’,STATUS= , UNKNOWN’) 

DO 80 I=1,NIDP1 
WRITE(12,110) I,D(I),FS(I) 

80 CONTINUE 

******************************************************************** 

* FORMATS 

* 

******************************************************************** 

100 FORMAT! IX, 15) 

110 FORMAT! IX, 15, 2X,F8.4,2X,F8.4) 

1001 FORMAT! 1X,I5,2X,F8.4,2X,F8.4,2X,E 10.4) 

1000 FORMAT! 1X,I5,2X,E 16. 10,2X,E16. 10) 

1010 F0RMAT(1X,I5,2X,I5,2X,F8.4) 

******** e********* ******************* ************************ 

END 

B. MAXIMUM LIKELIHOOD ESTIMATOR 

PROGRAM MLE 
INTEGER NIA 
PARAMETER (NLA = 6000) 

REAL FE(NIA),DFE(NIA),FUB(NLA),DFUB(NLA),GL 
CHARACTER* 1 Ql 

* 

PRINT *, ’Do you want the upper bound or complete data’ 
PRINT *, ’set, suggest the compete data set when the’ 

PRINT *, ’DATA is less the the number of filements’ 

110 PRINT *, TYPE "U" for upper and ”C" otherwise’ 

303 FORMAT(Al) 

READ (*,303) Ql 
IF(Ql .EQ. ’IT) THEN 

OPEN(UNIT= 1 1 ,FILE=’UBFECDF.OUT,STATUS=’OLD’) 
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ELSEIF (Q1 .EQ. ’C’) THEN 

OPEN(UNIT= 1 1 ,FILE=TECDF.OUT,STATUS=’OLD’) 

ELSE 

PRINT *, TRY TYPING CAPITAL LETTERS’ 

GOTO 110 
END IF 

* 

READU1.1010) M,N,EB 
1010 F0RMAT(1X,I5,2X,I5,2X,F8.4) 

DO 10 1=1, M 

READ(11,1020) J,DFUB(I),FUB(I) 

10 CONTINUE 

1020 FORMAT(lX,I5,2X,F8.4,2X,F8.4) 

************************************************************** 

* CALCULATE ALPHA AND BETA 

PRINT *, What is the GAGE LENGTH?’ 

READ *, GL 

* 

* REDUCE THE DATA TO ABOUT 50 POINTS 
888 PRINT *, TIow many points for the MLE ?’ 

READ *, MLE 
PRINT *, MLE 
DEL = 1.0/(MLE+1.0) 

J = 1 
I =1 

* find the value of F above and below the expected rank and average 
400 FE(I) = 1.0*I/(MLE+1.0) 

IF(I .GT. MLE) GOTO 450 

410 IF((FUB(J) .LE. FE(I)) .AND. (FUB(J+1) .GT. FE(I))) THEN 

DFE(I)= -(DFUB(J+1)-DFUB(J))*(FUB(J+1)-FE(I))/(FUB(J+1)-FUB(J)) 
+ + DFUB(J+1) 

1 = 1+1 

GOTO 400 

ELSEIF((J+1) .LE. (M)) THEN 

J=J+1 

GOTO 410 

ENDIF 

450 PRINT *, I,J,FUB(J),FUB(J+1),FE(I) 

IF(FUB(J+1) .LE. 0.0) THEN 
PRINT *, ’MLE to large try again’ 

GOTO 888 
ENDIF 

***** ** *** ****** ******** ************ ***************** ********* 

* estimate the initial alpha 

PRINT *, ’Estimate the intial alpha’ 

READ *, ALFA0 

************************************************************** 

* calculate 

ALFA = ALFA0 
KK-0 

320 DO 310 IT = 1,2 
SUM1=0.0 
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SUM2=0.0 

SUM3=0.0 

* Calculate the variables of alpha 

DO 300 1=1, MLE 
IF(DFE(I) .GT. 0.0) THEN 

SUM1 = SUMl + ((DFE(I))**ALFA * LOG(DFE(I))) 

SUM2 = SUM2 + ((DFE(I))**ALFA) 

SUM3 = SUM3 + LOG{DFE(I)) 

END IF 

300 CONTINUE 

IF(IT .EQ. 1) THEN 

ALFA1=1.0/(SUM1/SUM2 - (SUM3/MLE)) 

ALFA = ALFA1 
ELSEIFCIT .EQ. 2) THEN 
ALFA2=1.0/(SUM1/SUM2 - (SUM3/MLE)) 

ENDIF 

310 CONTINUE 

ALFAOP = ALFA0 -((ALFAO - ALFA1)**2)/(ALFA0 - 2*ALFA1 
+ + ALFA2) 

TEST = ABS(ALFAOP-ALFAO) 

IF(TEST .GT. .0001) THEN 
IF(KK .GT. 1000) GOTO 341 
KK = KK + 1 
ALFAO = ALFAOP 
PRINT *, ’ALFAO =’, ALFAO 
GOTO 320 

T?T QT? 

ALPHA = ALFAOP 

ENDIF 

SUM = 0.0 

DO 340 I = 1,MLE 

SUM = SUM + ( DFE( I ))** ALPHA 

340 CONTINUE 

341 BETA = (SUM/(MLE))**( 1/ALPHA) 

PRINT *, ’ALPHA = ’,ALPHA,’ BETA = ’,BETA/GL,’ mm/mm’ 
END 

C. WEIBULL UNEAR CDF 

PROGRAM LIN 
INTEGER NLA, AREA 
PARAMETER (NIA = 6000) 

REAL FE(NIA),DFE(NIA),FUB(NIA),DFUB(NIA),FS(NIA),D(NIA) 

* 

OPEN(UNIT=ll,FILE=’UBFECDF.OUT’,STATUS=’OLD’) 
READU1.1010) M,N,EB 
1010 F0RMAT(1X,I5,2X,I5,2X,F8.4) 

DO 10 1=1, M 

READ(1 1,1020) J,DFUB(I),FUB(I) 

10 CONTINUE 

1020 F0RMAT(1X,I5,2X,F8.4,2X,F8.4) 

************************************************************** 
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* REDUCE THE DATA TO ABOUT 50 POINTS, USING THE EXPECTED 
RANK 

PRINT *, HOW MANY POINTS TO FILTER USING THE EXPECTED 

RANK 

READ *, NUM 
PRINT *, NUM 
DEL = 1.0/(NUM+1.0) 

J = 1 
I =1 

* find the value of F above and below the expected rank and average 
400 FE(I) = 1.0*I/(NUM+1.0) 

IF(I .GT. NUM) GOTO 450 

410 IF((FUB(J) .LE. FE(I» .AND. (FUB(J+1) .GT. FE(I))) THEN 

DFE(I)= -(DFUB(J+1)-DFUB(J))*(FUB(J+1 )-FE(I ))/(FUB( J+ 1 )-FUB(J)) 

+ + DFUB(J+1) 

1 = 1+1 

GOTO 400 

ELSEIF((J+1) .LT. (M+l)) THEN 

J=J+1 

GOTO 410 

ENDIF 

450 PRINT *, I,J 

******************************************************************** 

* CALCULATE ALPHA AND BETA 

* 

* FIND F*=LOG{-LOG{l-F)) VS LOG{DELTA/BETA) 

DO 20 1=1, M 

IF(FUB(I) .LE. 0.0) FUB(I) = 0.0001 
FS(I) = LOG(-LOG{ 1-FUB(I ))) 

IF(DFUB(I) .LE. 0.0) DFUB(I) = 0.0001 
Dd) = LOG<DFUB(I)) 

20 CONTINUE 

DO 70 1=1, M 

510 IF((FS(I) .LE. 0.63212) .AND. (FSd+1) -GT. 0.63212)) THEN 
BETA = -(Dd+l)-Dd))*(FS(I+l)-0.63212)/(FS(I+l)-FS(I)) 

+ + Dd+1) 

GOTO 71 
ENDIF 

70 CONTINUE 

71 CONTINUE 
BETA = EXP(BETA) 

PRINT *, ’ BETA = \BETA 



DO 80 1=1, M 

IF(DFUB(I) .LE. 0.0) DFUB(I) = 0.0001 
D(I) = LOG{DFUB( I)/BETA) 

80 CONTINUE 

BETA = EB*BETA/N 

******************************************************************** 

OPENtUNIT^l.FILE^FSTAR.OU'r.STATUS^UNKNOWN’) 

DO 60 1=1, M 
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WRITE(11,1020) I,D(I),FS(I) 

CONTINUE 

END 
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APPENDIX F 



SIMULATION FORTRAN CODE 
PROGRAM SIMUL 



SIMUL IS A STATISTICAL SIMULATION ROUTINE USED 
FOR RELIABILITY 

ANALYSIS OF COMPOSITE FILEMENTS. 



UTILITIES FILES REQUIRED ARE: 

SEED - file containing the integer seed for random 
number generation 

INPUT - file containing the input parameters * 



OUTPUT FILES INCLUDE: 

EXPER.OUT - the experimental simulation of discrete 
data. 

SEED - outputs a new seed. 



INPUT FILE FOR BUNDLE2.FOR PROGRAM 
XXXXXXXXXXXXXXXXXXX12345678 
E (kg) .90 

ROWS = M 100 

COLS = N 100 

MGL (mm) 250.0 

ALPHA 4.0 

BETA 16.0 

MODEL NORMAL 

SKEW/ALPHA 3.5 

MAX SLACK (% ) 10.0 

SAMPLE RATE 20.0 

CROSS-HEAD SPEED 1.8 
NOISE +/- kg .001 

TOLD (mm) .00254 



dine 1 
dine 2 

:line 3 
dine 4 
:line 5 
dine 6 
dine 7 
:line 8 

dine 9 
dine 10 

dine 11 % OF MAX 
DISPLACEMENT 
:line 12 pts/s 
dine 13 mm/min 
:line 14 tolerance 
:line 15 tolerance of displ. 
output 



INPUT PARAMETERS ARE: 

1) THE LENGTH OF THE TEST BUNDLE (MGL). 

2) THE NUMBER OF FILEMENTS PER BUNDLE (N). 

3) THE NUMBER OF SEGMENTS PER FILEMENT (M). 

4) THE SEGMENT LENGTH IS DEFINED BY MGL=M*LS 



THIS PROGRAM CREATES: 



* 1) THE RANDOM NUMBERS 

* 2) CONVERTS TO WIEBULL BUNDLE 

* 3) FIND THE WEAKEST IN THE ENTIRE BUNDLE 

* FOR THE LOWER BOUND. 

* 4) FIND THE WEAKEST IN EACH FILEMENT 

* 5) ORDER THE WEAKEST TO STRONGEST FILEMENTS 

* 6) CHANGE THE STRENGTHS TO BREAKING LOADS AND 

* DISPLACEMENTS 

*************************************************************** 

* Begin Program: 

* 

*************************************************************** 



INTEGER ROWS, COLS, MAX, LIM.LIMl ,N,PTS, COUNT, LIM2 
PARAMETER! LIM=6000,LIM1=1000,LIM2=1001) 

REAL BUN(LIM1),RAND1(LIM1),RAND2(LIM1),DISPL(LIM1), 
+SLACK(LIM1),P(LIM1), ALPHA, BETA, MGL,BRKLDM, 
+DISPLM,MAXX,MAXY,MEAN,PU(LIM1 ), 
+PL(LIM1),SUM,MAXSLK,TEMP(LIM2),PS(LIM2),D(LIM), 
+DISPLA(LIMl),DEL,PMS,LOAD(LIM),NSE 
DIMENSION N0(LIM1),NB(LIM1),NDB(LIM1) 

DOUBLE PRECISION SEED 
CHARACTER* 1 MODEL*7,SONS,Q,Q2,Q3,Q4 
COMMON MAX 
COMMON /SML/M.N 

* 

******************************************************************** 



* FILE 11 = PARAM DATA FOR INPUT 

* FILE 12 = SEED DATA FOR SEEDING 

* 

******************************************************************** 



OPEN(UNIT=ll,FILE=’INPU'F,STATUS= , OLD’) 

OPEN(UNIT=12,FILE=’SEED’,STATUS=’OLD’) 

* 

***^c*4:********4;****4:*4:*4:4:**4:******4;4:4:*4:*4;4c******************** 

* READ IN THE INPUT DATA 

* 

************************************************************** 



READ! 11, 1530) 

READ! 11, 1500) E 
PRINT *, ’MODULUS^, E 
READ! 11,1520) M 
PRINT *, *M = ’,M 
ROWS=M 
READ(11,1520) N 
COLS=N 

PRINT *, *N = ’,N 

READ(11,1510) MGL 

PRINT*, ’Mean Gage Length = ’,MGL 
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MAX=M*N 

READd 1,1510) ALPHA 
PRINT *, ’ALPHA = ALPHA 
READ(11,1510) BETA 
PRINT *, 'BETA =’,BETA 

* 

****** 

* Determine which model to use for slack. 
****** 



READ(11,1515) MODEL 

PRINT *, ’MODEL=’,MODEL,’ APPROXIMATION 
IF (MODEL .EQ. NORMAL’) THEN 
READ(11,1500) STDDEV 

PRINT *, ’ALPHA = ’, STDDEV 
REATK1L1510) MEAN 

PRINT *, ’Maximum slack as percent of expected’ 

PRINT *, ’maximum displacement (3% gage length). =’, ME AN 
MEAN = (MEAN/1 00)*(0.03)*MGIV2 
PRINT *, ’Mean=\MEAN,’mm’ 

PRINT *, ’Maximum slack =’,MEAN*2,’mm’ 

ELSE 

READdl.lSOO) RANGE1 
READdl,1500) RANGE2 
RANGE1=(RANGE1/100)*(.03)*MGL 
RANGE2=(RANGE2/100)*(.03)*MGL 
PRINT *, ’RANGE 1 =’, RANGE 1,’ mm’ 

PRINT *, ’RANGE 2 =’,RANGE2,’ mm’ 

ENDIF 

* 

****** 

* Read in the sampling rate and cross-head speed. 
****** 

***********************5* e******************************************* 

* Ensure the displacement interval governed by the cross-head and 

* sampling rate is greater than the tolerance level of the output 

* displacement data. 

* 

************************* ******** ************************ ********** 



READdl,1510) SR 
PRINT *, ’Sampling rate =’,SR 
READ(11,1510) XSP 
PRINT *, ’Cross-head speed =’,XSP 

* 

****** 
* How much noise. 

****** 



READ(11,1510) NSE 
PRINT *, Ttange of noise =’,NSE 
READ(11,1510) TOLD 
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PRINT *, Tolerance of output =’,TOLD 
********* 

* Calculate the number of total points recorded in the * 

* experiment. * 

********* 

* 

444 RR = MGL*(0.03) 

PTS = RR*SR/(XSP/60) 

DELTAD = RR/PTS 
IF(TOLD .GT. DELTAD) THEN 

PRINT *, The displacement interval is less than the output’ 
PRINT *, ’tolerance, please input a new cross-head and/or’ 
PRINT *, ’sampling rate.’ 

PRINT *, ’Current Cross-head = ’ r XSP 
PRINT *, ’Input new XSP’ 

READ *, XSP 

PRINT *, ’Current sampling rate = ’,SR 
PRINT *, ’Input new sampling rate’ 

READ *, SR 
GOTO 444 
ENDIF 

* 

****** 

* Determine the type of simulation. 
****** 



PRINT *, ’SLACK OR NO SLACK ? (Y/N)’ 

READ (*,1600) SONS 

PRINT *, ’NOISE OR NO NOISE ? (Y/N)’ 

READ (*,1600) Q2 

PRINT *, ’Change the SEED ?’ 

READ (*,1600) Q 

PRINT *, ’COMPLAINCE (Y/N)?’ 

READ (*,1600) Q3 

IF((Q3 .EQ. T) .OR. (Q3 .EQ. ’y’)) THEN 

PRINT *, ’WHAT IS THE COMPLIANCE NOISE LEVEL?’ 

READ *, CTOL 

ENDIF 

PRINT *, T)o you want English units (lbs.inches) on output? 
+(Y/N)’ 

PRINT *, ’If N the output is in SI units (kg, mm)’ 

READ (*,1600) Q4 

CLOSE(ll) 
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******************************************************************** 



* BEGIN THE ALGORITHM 

^^^^^3tc^^^^^^3^3|c3fc3|c3fc3|ca|ca|c3|c3|c3|c3|c3|c9|c9|caf:a|c9«ca|catcaf:9|c3|ca(ca|c9«ca4ca|ca|c3|c9|ca|c9|c9|ca«caf:a|caf:****************** 

** 

********* 

* Read the seed, (seed for the uniform random number generator). 

* and generate the random numbers 
********* 

* 

READ( 12,1000) SEED 
REWINIX12) 

DO 10 1=1, N 

CALL RAND(SEED,M,BUN) 

CALL SMALL(R0WS,BUN,RAND1(I)) 

10 CONTINUE 



********* 

* Check if slack is requested in simulation. 
********* 



IF (SONS .EQ. TO THEN 
CALL RAND(SEED,N,RAND2) 
ENDIF 



********* 
* Initialize the order arrays. 

** ******* 



DO 25 I=1N 
SLACK(I)=0.0 
NO(I)=I 
NB(I)=I 
NDB(I)=I 

25 CONTINUE 



********* 

* Check if seed is to be changed, if so write the new seed. 
********* 



* 

IF ((Q .EQ. TO .OR. (Q .EQ. ’yO) THEN 
WRITE( 12,1000) SEED 
REWIND (12) 

CLOSE(12) 

ENDIF 

* 

********************************* *********************************** 

* Generate the load and slack distributions. 

* 

******************************************************************** 

* 
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CALL WEIBUL(ALPHA,BETA,N,RAND1,P) 

IF (SONS .EQ. T* .OR. SONS .EQ. V) THEN 
IF (MODEL .EQ. ’NORMAL’) THEN 

CALL NORMAL(STDDEV,MEAN,N,RAND2, SLACK, MAXSLK) 
ELSE 

CALL UN IFM(RANGE1,RANGE2,N,RAND2, SLACK, MAXSLK) 
ENDIF 
SLACK(1)=0.0 
ENDIF 

* 

******************************************************************** 

* Find the displacement for each failure load and slack. 

* 

******************************************************************** 



DO 20 1=1, N 

DISPLd) = (P(I)/E)*(SLACK(I) + MGL) 
20 CONTINUE 
K=N 

****** 

* Order the "DISPL" array, and carry "No” 
****** 



CALL SHELUN, DISPL, NO) 

* 

****** 

* reorder the remaining arrays, the others are 

* keyed to "DISPL" 

****** 



CALL SWITCH(N,NO,RANDl) 

CALL SWITCH(N r NO,RAND2) 

CALL SWITCH(N,NO, SLACK) 

CALL SWITCH(N,NO,P) 

* 

******************************************************************** 

* Calculate the apparent displacement. 
******************************************************************** 



DO 40 1=1 ,N 

DISPLA(I) = DISPLd) + SLACK(I) 

40 CONTINUE 
* 

******* 

* Reorder the arrays keyed to "DISPLA" 
******* 



CALL SHELL(N, DISPLA ,NB) 
CALL SWITCH(N,NB,RAND1) 
CALL SWITCH(N r NB,RAND2) 
CALL SWITCH(N ,NB , SLACK) 
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